% Function to compute the equilibrium of the model, considering selective
% assignment of students
function equilibrium = EquilibriumFunction_SelectiveAssignment_Counterfactual(parameter_vector_all,Pop_data,pi_E,EE_set,HousingCostOption,NE,Nxr,Nxn,paramsR,paramsF,paramsC,xr,xn,Xn,Xr,XrXr,XnXn,frn_e,index_educ)

alphaS=parameter_vector_all(1);
cSbar=parameter_vector_all(2);
cAbar=parameter_vector_all(3);
beta=parameter_vector_all(4);
gamE=parameter_vector_all(5);
etaE=parameter_vector_all(6);  
tau_F=parameter_vector_all(7);
xiH0_F=parameter_vector_all(8);
xiH0_C=parameter_vector_all(9);
XA=parameter_vector_all(10);
XS=parameter_vector_all(11);
XM=parameter_vector_all(12);
xiH1_F=parameter_vector_all(13);
xiH1_C=parameter_vector_all(14);
varphi=parameter_vector_all(15);
alphaA=parameter_vector_all(16);
tau_F_income=parameter_vector_all(17);
gam=parameter_vector_all(18);
tau_R=parameter_vector_all(19); 
tau_C=parameter_vector_all(20);
tau_1_F=parameter_vector_all(21);
alphaM = 1 - alphaA - alphaS;
etaE_R =  etaE;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Big Loop: Computing the Equilibrium Locations and Equilibrium Prices, Quantities
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Parameters of the Big Loop (over locations)    
        dist    =   10;
        J       =   1;
        Jmax    =   500;
        rlx     =   0.1;
        tol     =   0.001;
        
%% Some (initial or guessed) distribution of the population
        Pop_R = Pop_data(1);
        Pop_F = Pop_data(2);
        Pop_C = 1 - Pop_R - Pop_F;
        PoP   = 1;

%% Location Choices: An Initial guessed value to start the iterations
        q_R0  = Pop_R*ones(Nxn,Nxr,NE);
        q_F0  = Pop_F*ones(Nxn,Nxr,NE);
        q_C0  = Pop_C*ones(Nxn,Nxr,NE);
        q_U0  = q_F0 + q_C0;

        q_R  = q_R0;
        q_F  = q_F0;
        q_C  = q_C0;
        q_U  = q_F0 + q_C0;

    
%% ---Speeding Up the Iterations, as the loop closes in the convergence
%  ---use as the initial values in the fminsearch functions.
        p_SR0 = alphaS;
        p_SU0 = alphaS/alphaA;    
        p_M0  = alphaM/alphaA;    

        p_SR  = p_SR0;
        p_SU  = p_SU0;
        p_M  = p_M0;
 
%% Constructing the implied skill distribution across regions, given the country's education distribution, piE, the location choices, q_l, and the conditional skill distribution, frn_e(:,:,e) 
        while J<Jmax;      
        % Comment: I added the term 1/PoP to make sure that the population adds to one; numerical integrals only get it close to 1. I only cor it once, at the first iteration.     
        if NE==5,
            frn_R   = 1/PoP*( pi_E(1)*frn_e(:,:,1).*q_R(:,:,1) + pi_E(2)*frn_e(:,:,2).*q_R(:,:,2) + pi_E(3)*frn_e(:,:,3).*q_R(:,:,3) + pi_E(4)*frn_e(:,:,4).*q_R(:,:,4) + pi_E(5)*frn_e(:,:,5).*q_R(:,:,5) );
            frn_F   = 1/PoP*( pi_E(1)*frn_e(:,:,1).*q_F(:,:,1) + pi_E(2)*frn_e(:,:,2).*q_F(:,:,2) + pi_E(3)*frn_e(:,:,3).*q_F(:,:,3) + pi_E(4)*frn_e(:,:,4).*q_F(:,:,4) + pi_E(5)*frn_e(:,:,5).*q_F(:,:,5) );
            frn_C   = 1/PoP*( pi_E(1)*frn_e(:,:,1).*q_C(:,:,1) + pi_E(2)*frn_e(:,:,2).*q_C(:,:,2) + pi_E(3)*frn_e(:,:,3).*q_C(:,:,3) + pi_E(4)*frn_e(:,:,4).*q_C(:,:,4) + pi_E(5)*frn_e(:,:,5).*q_C(:,:,5) );
        elseif NE==4,
            frn_R   = 1/PoP*( pi_E(1)*frn_e(:,:,1).*q_R(:,:,1) + pi_E(2)*frn_e(:,:,2).*q_R(:,:,2) + pi_E(3)*frn_e(:,:,3).*q_R(:,:,3) + pi_E(4)*frn_e(:,:,4).*q_R(:,:,4));
            frn_F   = 1/PoP*( pi_E(1)*frn_e(:,:,1).*q_F(:,:,1) + pi_E(2)*frn_e(:,:,2).*q_F(:,:,2) + pi_E(3)*frn_e(:,:,3).*q_F(:,:,3) + pi_E(4)*frn_e(:,:,4).*q_F(:,:,4));
            frn_C   = 1/PoP*( pi_E(1)*frn_e(:,:,1).*q_C(:,:,1) + pi_E(2)*frn_e(:,:,2).*q_C(:,:,2) + pi_E(3)*frn_e(:,:,3).*q_C(:,:,3) + pi_E(4)*frn_e(:,:,4).*q_C(:,:,4));
        end

        for ie=1:NE,
            Pi_E_R(ie)=pi_E(ie)*trapz(xr,trapz( xn,frn_e(:,:,ie).*q_R(:,:,ie) ) );
            Pi_E_F(ie)=pi_E(ie)*trapz(xr,trapz( xn,frn_e(:,:,ie).*q_F(:,:,ie) ) );
            Pi_E_C(ie)=pi_E(ie)*trapz(xr,trapz( xn,frn_e(:,:,ie).*q_C(:,:,ie) ) );
        end

    % Normalizaing to one, i.e. cross-section education distribution in each region
        pi_E_R =Pi_E_R./sum(Pi_E_R);
        pi_E_F =Pi_E_F./sum(Pi_E_F);
        pi_E_C =Pi_E_C./sum(Pi_E_C);

%% Education distribution
        educ_dist_R = Pop_R*pi_E_R;
        educ_dist_F = Pop_F*pi_E_F;
        educ_dist_C = Pop_C*pi_E_C;
      
%% Average Education in each Region
        index_educ_op = 5-index_educ;
        EE_R = pi_E_R*EE_set;
        EE_F = (educ_dist_F(1:index_educ_op)*EE_set(1:index_educ_op))/sum(educ_dist_F(1:index_educ_op));
        EE_C = (educ_dist_F((index_educ_op+1):5)*EE_set((index_educ_op+1):5)+educ_dist_C*EE_set)/(sum(educ_dist_F((index_educ_op+1):5))+sum(educ_dist_C));
    

%% Constructing the Conditional Probabilities      
        for ie0=1:NE,
            % Conditional Transition Probabilities
                Pr_ee_R(ie0,:) = TransProb(ie0,EE_R,paramsR)';
                Pr_ee_F(ie0,:) = TransProb(ie0,EE_F,paramsF)';
                Pr_ee_C(ie0,:) = TransProb(ie0,EE_C,paramsC)';
            % Conditional Expectation of Children's Education, for each household, by Region
                CEE_R(ie0) = Pr_ee_R(ie0,:)*EE_set;
                CEE_F(ie0) = Pr_ee_F(ie0,:)*EE_set;
                CEE_C(ie0) = Pr_ee_C(ie0,:)*EE_set;
        end

        if index_educ==1
            Pr_ee_F(5,:) = TransProb(5,EE_C,paramsF)';
            CEE_F(5) = Pr_ee_F(5,:)*EE_set;
        end

        if index_educ==2
            Pr_ee_F(5,:) = TransProb(5,EE_C,paramsF)';
            CEE_F(5) = Pr_ee_F(5,:)*EE_set;
            Pr_ee_F(4,:) = TransProb(4,EE_C,paramsF)';
            CEE_F(4) = Pr_ee_F(4,:)*EE_set;
        end

        if index_educ==3
            Pr_ee_F(5,:) = TransProb(5,EE_C,paramsF)';
            CEE_F(5) = Pr_ee_F(5,:)*EE_set;
            Pr_ee_F(4,:) = TransProb(4,EE_C,paramsF)';
            CEE_F(4) = Pr_ee_F(4,:)*EE_set;
            Pr_ee_F(3,:) = TransProb(3,EE_C,paramsF)';
            CEE_F(3) = Pr_ee_F(3,:)*EE_set;
        end
        if index_educ==4
            Pr_ee_F(5,:) = TransProb(5,EE_C,paramsF)';
            CEE_F(5) = Pr_ee_F(5,:)*EE_set;
            Pr_ee_F(4,:) = TransProb(4,EE_C,paramsF)';
            CEE_F(4) = Pr_ee_F(4,:)*EE_set;
            Pr_ee_F(3,:) = TransProb(3,EE_C,paramsF)';
            CEE_F(3) = Pr_ee_F(3,:)*EE_set;
            Pr_ee_F(2,:) = TransProb(2,EE_C,paramsF)';
            CEE_F(2) = Pr_ee_F(2,:)*EE_set;
        end
  

%% Some other auxiliary distributions   
        f_XrXn_R      = frn_R./XrXr./XnXn;
        Xr_f_XrXn_R   = frn_R./XnXn;
        Xn_f_XrXn_R   = frn_R./XrXr;

        f_XrXn_F      = frn_F./XrXr./XnXn;
        Xr_f_XrXn_F   = frn_F./XnXn;
        Xn_f_XrXn_F   = frn_F./XrXr;

        f_XrXn_C      = frn_C./XrXr./XnXn;
        Xr_f_XrXn_C   = frn_C./XnXn;
        Xn_f_XrXn_C   = frn_C./XrXr;
    
        f_XrXn_U              =    f_XrXn_C +    f_XrXn_F;
        Xr_f_XrXn_U           = Xr_f_XrXn_C + Xr_f_XrXn_F;
        Xn_f_XrXn_U           = Xn_f_XrXn_C + Xn_f_XrXn_F;
    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Equilibrium, conditional on Location Decisions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    wu=XA; % pA=1, numeraire

%% Rural Regions
% Setting up the function that prices services and allocates labor in Rural Areas
% Populations
    pop_R = trapz( Xr,trapz(Xn,   f_XrXn_R ));  % employment in routine tasks, rural areas
% calling the function to find the equilibrium in Rural Regions
    p_SR0=p_SR; % After first using an arbitrary initial value, just use the previous equilibrium values 

    fun0=@(rural_p) RuralEquilibriumFunction(rural_p,pop_R, XS,XA,alphaA,alphaM,alphaS,cAbar,cSbar, Xn,Xr,XrXr,XnXn,f_XrXn_R);    
        options = optimset('TolFun',1e-8);
        options = optimset('MaxIter',2000);
        options = optimset('MaxFunEvals',5000);
    [p_SR,fval_rural,exitflag_Rural,output] = fminsearch(fun0,p_SR0,options);

%Just Making Sure it is positive:
    p_SR=max(p_SR,0.001);
      
%% Equilibrium in Rural Areas: Implied wages and labor market earnings
    wrR                     =    p_SR*XS;
    ybar_R                  =    p_SR*cSbar*(1-alphaS)/alphaS+cAbar;       % income threshold for the consumption of Services (presumes zero housing costs in rural areas)
% Labor Markets
    yR                      =    max(wu,wrR*XrXr);
    IOCC_r_R                =    zeros(size(XrXr));  % indicator, occupation choices
    IOCC_r_R(wu<wrR*XrXr)   =    1;
    ICS_R                   =    zeros(size(XrXr));  % Indicator, positive 
    ICS_R(yR>=ybar_R)       =    1;
    LrR                     =    trapz( Xr,trapz(Xn, Xr_f_XrXn_R .*IOCC_r_R));  % effective labor in routine tasks, rural areas
    ErR                     =    trapz( Xr,trapz(Xn,    f_XrXn_R .*IOCC_r_R));  % employment in routine tasks, rural areas
    QA                      =    XA*(pop_R-ErR);
    QS_R                    =    XS*LrR;
% Consumption of Services
    yTop_R                  =    trapz( Xr,trapz(Xn, yR.*f_XrXn_R.*ICS_R));
    fTop_R                  =    trapz( Xr,trapz(Xn,     f_XrXn_R.*ICS_R));
% Equilibrium Pricing Condition
    p_SR1                  =    alphaS*(yTop_R-cAbar*fTop_R)/(QS_R+(1-alphaS)*cSbar*fTop_R);
    QA_eq           =   QA;
    PSR_eq          =   p_SR;


%% Urban Regions   
%% Finding the No-Services equilibrium conditions in Urban regions
    rw_rn0 =  varphi/(1-varphi);
    funM=@(wr_wn) M_OnlyFunction(wr_wn,varphi,Xr,Xn,XrXr,XnXn, f_XrXn_U, Xr_f_XrXn_U, Xn_f_XrXn_U );
        options = optimset('TolFun',0.00000000000001);
        options = optimset('MaxIter',1000);
        options = optimset('MaxFunEvals',5000);
    [rw_rn,fval_Monly,exitflag_M_Only,output] = fminsearch(funM,rw_rn0,options);

% just making sure it is positive....
    rw_rn=max(0.001,rw_rn);
    
% M only allocations
    IOCC_r_U_M_only                        =  zeros(size(XrXr));
    IOCC_n_U_M_only                        =  zeros(size(XrXr));
    IOCC_r_U_M_only(rw_rn*XrXr >  XnXn)    =  1;
    IOCC_n_U_M_only(rw_rn*XrXr <  XnXn)    =  1;
    L_r_M_only                      =  trapz( Xr,trapz(Xn,  Xr_f_XrXn_U.*IOCC_r_U_M_only));
    Emp_r_U_M_only                  =  trapz( Xr,trapz(Xn,     f_XrXn_U.*IOCC_r_U_M_only));
    L_n_M_only                      =  trapz( Xr,trapz(Xn,  Xn_f_XrXn_U.*IOCC_n_U_M_only));
    Emp_n_U_M_only                  =  trapz( Xr,trapz(Xn,     f_XrXn_U.*IOCC_n_U_M_only));
    wrn_M_only                      =  rw_rn;
    Q_M_only                        =  XM*(L_r_M_only)^varphi*(L_n_M_only)^(1-varphi);
    pSM_M_only                      =  (XM/XS)*varphi*(L_n_M_only/L_r_M_only)^(1-varphi);


%% Setting up the function that solves for Urban Equilibrium (prices and allocations)
% Populations
    pop_C                   =    trapz( Xr,trapz(Xn,   f_XrXn_C ));  % employment in routine tasks, rural areas
    pop_F                   =    trapz( Xr,trapz(Xn,   f_XrXn_F ));  % employment in routine tasks, rural areas
    pop_U                   =    pop_C + pop_F;
    pop = [pop_R; pop_F; pop_C];
    pop = pop./sum(pop);
    pi_E_U = (pi_E_F.*pop_F+pi_E_C.*pop_C)./pop_U;

% Hosing Costs: in units of M goods required to construct a house
    xiH_R    =   0; 
    if HousingCostOption==1,
        xiH_F    =    xiH0_F*(1+pop_F)^xiH1_F;
        xiH_C    =    xiH0_C*(1+pop_C)^xiH1_C;
    elseif HousingCostOption==2,
        xiH_F    =   xiH0_F*(pop_F)^xiH1_F;
        xiH_C    =   xiH2_C+xiH0_C*(pop_C)^xiH1_C;
    elseif HousingCostOption==3,
        xiH_F    =   xiH0_F*exp(xiH1_F * pop_F);
        xiH_C    =   xiH2_C+xiH0_C*exp(xiH1_C * pop_C);
                
    end
        
        
%% Calling the Function for Urban Equilibrium 
% Using equilibrium prices from previous location's iteration as initial condition 
    p_SU0 = p_SU;    
    p_M0  = p_M;
    pp0=[p_SU0,p_M0]; 
% Calling the function and setting its options 
    fun1=@(urb_p) UrbanEquilibriumFunction(tau_F_income,urb_p,pop_R,pop_R,xiH_F,xiH_C,XS,XM,varphi,QA_eq,alphaA,alphaM,alphaS,cAbar,cSbar,Xn,Xr,XrXr,XnXn,f_XrXn_C,f_XrXn_F,L_r_M_only,L_n_M_only,Q_M_only,pSM_M_only);
        options = optimset('TolFun',1e-9);
        options = optimset('MaxIter',2000);
        options = optimset('MaxFunEvals',5000);
        options = optimset('TolX',1e-10);
    [urban_prices,fval_urb1,exitflag_Urban,output] = fminsearch(fun1,pp0,options);

%% Recovering the Equilibrium    
    p_SU    =  max(0.001,  urban_prices(1));
    p_M     =   max(0.001, urban_prices(2));
  
%% Constructing measures of skills distributions in Urban Regions (U= C + F)    
   f_XrXn_U   =  f_XrXn_C +f_XrXn_F;  
    Xn_f_XrXn_C   =  XnXn.*f_XrXn_C;
    Xn_f_XrXn_F   =  XnXn.*f_XrXn_F;
    Xn_f_XrXn_U   =  XnXn.*f_XrXn_U;
    Xr_f_XrXn_C   =  XrXr.*f_XrXn_C;
    Xr_f_XrXn_F   =  XrXr.*f_XrXn_F;
    Xr_f_XrXn_U   =  XrXr.*f_XrXn_U;

    factor_wn  =  varphi^( varphi/(1-varphi) )*(1-varphi)*(XM/XS^varphi)^(1/(1-varphi));

    if p_SU/p_M < pSM_M_only,
        wrU     =       varphi*p_M*XM*(L_n_M_only/L_r_M_only)^(1-varphi);
        wnU     =   (1-varphi)*p_M*XM*(L_n_M_only/L_r_M_only)^(-varphi);
    else
        wrU                             =    p_SU*XS;
        wnU                             =    factor_wn*( p_M /(p_SU)^varphi )^(1/(1-varphi));  % these are the wages implied by the case when both S and M are produced at U
    end

%% Individual Earnings
    yU                              =    max(wrU*XrXr, wnU*XnXn);   

%% Consumption of Services
    ybar_F                          =    p_SU*cSbar*(1-alphaS)/alphaS + cAbar + p_M*xiH_F;        % income threshold for the consumption of Services, households in slums
    ybar_C                          =    p_SU*cSbar*(1-alphaS)/alphaS + cAbar + p_M*xiH_C;        % income threshold for the consumption of Services, households in cities
    ICS_F                           =    zeros(size(XrXr));  % Indicator, positive consumption, households in slums
    ICS_C                           =    zeros(size(XrXr));  % Indicator, positive consumption, households in cities
    ICS_F(yU*(1-tau_F_income)>=ybar_F)               =    1;
    ICS_C(yU>=ybar_C)               =    1;
% Aggregation: Consumption of Services
    yTop_F                  =    trapz( Xr,trapz(Xn, yU.*f_XrXn_F.*ICS_F));    % Slums
    fTop_F                  =    trapz( Xr,trapz(Xn,     f_XrXn_F.*ICS_F));
    yTop_C                  =    trapz( Xr,trapz(Xn, yU.*f_XrXn_C.*ICS_C));    % Cities
    fTop_C                  =    trapz( Xr,trapz(Xn,     f_XrXn_C.*ICS_C));
 
% Aggregation: Supplies of Skills and Employment    
    IOCC_r_U                        =    zeros(size(XrXr));  % indicator, occupation choices
    IOCC_r_U(wnU*XnXn<wrU*XrXr)     =    1;
    LrU                             =    trapz( Xr,trapz(Xn, Xr_f_XrXn_U .*IOCC_r_U));      %  effective labor in routine tasks, urban areas
    ErU                             =    trapz( Xr,trapz(Xn,    f_XrXn_U .*IOCC_r_U));      %  employment in routine tasks, urban areas
    LnU                             =    trapz( Xr,trapz(Xn, Xn_f_XrXn_U .*(1-IOCC_r_U)));  %  effective labor in non-routine tasks, urban areas
    EnU                             =    trapz( Xr,trapz(Xn,    f_XrXn_U .*(1-IOCC_r_U)));  %  employment in non-routine tasks, urban areas
    if p_SU/p_M<pSM_M_only,
        QM                              =   Q_M_only;
        QSU                             =   0; 
    else
        LrM                                       =   (varphi*XM*p_M/XS/p_SU)^(1/(1-varphi))*LnU;
        QM  =   XM*(LrM)^varphi*(LnU)^(1-varphi);
        QSU =   XS*(LrU-LrM);          
    end

    p_SU1  =    alphaS*( yTop_F*(1-tau_F_income)+yTop_C - fTop_C*(cAbar+p_M*xiH_C) - fTop_F*( cAbar+p_M*xiH_F ) )/( QSU +(1-alphaS)*cSbar*(fTop_F+fTop_C) );
    p_M1  =    alphaM/alphaA*( QA_eq - cAbar )/( QM -xiH_C*pop_C-xiH_F*pop_F );
   

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
%%    Attained Utilities, Average Education Levels per Region, Location Choices
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%% Consumption
    %% Indirect utilities: consumption for each region and skills levels
        Factor_Below  = (cSbar)^alphaS*(alphaA)^alphaA*(alphaM)^alphaM./(alphaA+alphaM)^(alphaA+alphaM);
        Factor_Above  = (alphaS)^alphaS*(alphaA)^alphaA*(alphaM)^alphaM;

    %% Indirect Utilities of Consumption: Baseline: low value assuming non-feasibility
        IndU_Cons_R=-9999*ones(Nxn,Nxr);
        IndU_Cons_F=-9999*ones(Nxn,Nxr);
        IndU_Cons_C=-9999*ones(Nxn,Nxr);

    %% Minimum Income to Subsist in each Region
        y_subs_R = cAbar;
        y_subs_F = cAbar + xiH_F*p_M;
        y_subs_C = cAbar + xiH_C*p_M;

    if gam==1,
    %% Utility in Each Region, assuming no consumption of market services
        poor_IndU_Cons_R = log( (1-tau_R)*Factor_Below*(yR-y_subs_R)./(p_M)^alphaM );
        poor_IndU_Cons_F = log( (1-tau_F)*(Factor_Below^(1-tau_1_F))*((yU*(1-tau_F_income)-y_subs_F).^(1-tau_1_F))./(p_M)^(alphaM*(1-tau_1_F)) );
        poor_IndU_Cons_C = log( (1-tau_C)*Factor_Below*(yU-y_subs_C)./(p_M)^alphaM );

    %% Utility in Each Region, assuming positive consumption of market services, {ICS_R, ICS_F, ICS_C=1}
        rich_IndU_Cons_R = log( (1-tau_R)*Factor_Above*(yR + p_SR*cSbar - y_subs_R)./( (p_M)^alphaM*(p_SR)^alphaS ) );
        rich_IndU_Cons_F = log( (1-tau_F)*(Factor_Above^(1-tau_1_F))*((yU*(1-tau_F_income) + p_SU*cSbar - y_subs_F).^(1-tau_1_F))./( ((p_M)^alphaM*(p_SU)^alphaS )^(1-tau_1_F)) );
        rich_IndU_Cons_C = log( (1-tau_C)*Factor_Above*(yU + p_SU*cSbar - y_subs_C)./( (p_M)^alphaM*(p_SU)^alphaS ) );    
    else
    %% Utility in Each Region, assuming no consumption of market services
        poor_IndU_Cons_R = 1/(1-gam)*( (1-tau_R)*Factor_Below*(yR-y_subs_R)./(p_M)^alphaM ).^(1-gam);
        poor_IndU_Cons_F = 1/(1-gam)*( (1-tau_F)*(Factor_Below^(1-tau_1_F))*((yU*(1-tau_F_income)-y_subs_F).^(1-tau_1_F))./((p_M)^alphaM^(1-tau_1_F)) ).^(1-gam);
        poor_IndU_Cons_C = 1/(1-gam)*( (1-tau_C)*Factor_Below*(yU-y_subs_C)./(p_M)^alphaM ).^(1-gam);

    %% Utility in Each Region, assuming positive consumption of market services, {ICS_R, ICS_F, ICS_C=1}
        rich_IndU_Cons_R = 1/(1-gam)*( (1-tau_R)*Factor_Above*(yR + p_SR*cSbar - y_subs_R)./( (p_M)^alphaM*(p_SR)^alphaS ) ).^(1-gam);
        rich_IndU_Cons_F = 1/(1-gam)*( (1-tau_F)*(Factor_Above^(1-tau_1_F))*((yU*(1-tau_F_income) + p_SU*cSbar - y_subs_F).^(1-tau_1_F))./( ((p_M)^alphaM*(p_SU)^alphaS )^(1-tau_1_F)) ).^(1-gam);
        rich_IndU_Cons_C = 1/(1-gam)*( (1-tau_C)*Factor_Above*(yU + p_SU*cSbar - y_subs_C)./( (p_M)^alphaM*(p_SU)^alphaS ) ).^(1-gam);
    end

    %% Attained Utilities: Consumption 
        %% First layer: Can the household subsist?   
            IndU_Cons_R(yR>y_subs_R) = poor_IndU_Cons_R(yR>y_subs_R);
            IndU_Cons_F(yU*(1-tau_F_income)>y_subs_F) = poor_IndU_Cons_F(yU*(1-tau_F_income)>y_subs_F);
            IndU_Cons_C(yU>y_subs_C) = poor_IndU_Cons_C(yU>y_subs_C);

        %% Second layer: Would the household buy Services in the market?
            IndU_Cons_R(ICS_R==1) = rich_IndU_Cons_R(ICS_R==1);
            IndU_Cons_F(ICS_F==1) = rich_IndU_Cons_F(ICS_F==1);
            IndU_Cons_C(ICS_C==1) = rich_IndU_Cons_C(ICS_C==1);

%% Children's Expected Education (CEE)
    %% Utility for CEE for each region and education levels            
        if gamE==1,
            for ie=1:NE,
                WCEE_R(ie) = log(CEE_R(ie));
                WCEE_F(ie) = log((CEE_F(ie)));           
                WCEE_C(ie) = log(CEE_C(ie));
            end
        else
            for ie=1:NE,
                WCEE_R(ie) = 1/(1-gamE)*(CEE_R(ie))^(1-gamE);
                WCEE_F(ie) = 1/(1-gamE)*((CEE_F(ie)))^(1-gamE);
                WCEE_C(ie) = 1/(1-gamE)*(CEE_C(ie))^(1-gamE);
            end
        end
    
%% Attained Utilities: U[Consumption] + betaE*CEE_l   
        for ie=1:NE,
            V_R(:,:,ie) = (1-beta)*IndU_Cons_R + beta*WCEE_R(ie);
            V_F(:,:,ie) = (1-beta)*IndU_Cons_F + beta*WCEE_F(ie);
            V_C(:,:,ie) = (1-beta)*IndU_Cons_C + beta*WCEE_C(ie);
        end


%% Equilibrium Location Decisions
        q_R1 = exp(etaE_R*V_R)./( exp(etaE_R*V_R) + exp(etaE*V_F) + exp(etaE*V_C));
        q_F1 = exp(etaE*V_F)./( exp(etaE_R*V_R) + exp(etaE*V_F) + exp(etaE*V_C));
        q_C1 = exp(etaE*V_C)./( exp(etaE_R*V_R) + exp(etaE*V_F) + exp(etaE*V_C));

        dist = max(max(max(abs(q_C-q_C1))))+ max(max(max(abs(q_F-q_F1))))+ max(max(max(abs(q_R-q_R1))));
    
        q_R= ( rlx*q_R1+(1-rlx)*q_R );
        q_F= ( rlx*q_F1+(1-rlx)*q_F );
        q_C= ( rlx*q_C1+(1-rlx)*q_C );
    
        pi_E_L = [pi_E_R; pi_E_F; pi_E_C]';
        pi_E_calc = pi_E_L*pop;
        pi_E_calc = pi_E_calc./sum(pi_E_calc);

        if dist<=tol,
            J=Jmax;
        else
            J=J+1;    
        end

end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Reporting!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Total Output, 
    GDP    =  QA+p_M*QM+p_SR*QS_R+p_SU*QSU;
    VA_A   = 100*QA/GDP;
    VA_SR  = 100*p_SR*QS_R/GDP;
    VA_SU  = 100*p_SU*QSU/GDP;
    VA_M   = 100*p_M*QM/GDP;
    
    GNI    =  wu*(pop_R-ErR)+wrR*LrR+wrU*LrU+wnU*LnU;
    
    PYS_u  = 100*wu*(pop_R-ErR)/GNI; 
    PYS_rR  = 100*wrR*LrR/GNI; 
    PYS_rU  = 100*wrU*LrU/GNI; 
    PYS_nU  = 100*wnU*LnU/GNI; 

format short


%% Results:

% Labor distribution:
     l_A = (pop_R-ErR);
     l_SR = ErR;
     if QSU==0
         l_SU=0;
     else
        l_SU = (1-LrM/LrU)*ErU;
     end
     l_S = l_SR + l_SU;
     if QSU==0
         l_M=ErU + EnU;
     else
        l_M = (LrM/LrU)*ErU + EnU;
     end
     sum_l = l_A+l_M+l_S;
     l_A = l_A/sum_l;
     l_M = l_M/sum_l;
     l_S = l_S/sum_l;
     
% Average education in urban areas:
     EE_U = (EE_F*(pop_F/(pop_F+pop_C))+EE_C*(pop_C/(pop_F+pop_C)));
     
% Average education and education distribution in the next period:
     EE_R_future = pi_E_R*CEE_R';
     EE_F_future = pi_E_F*CEE_F';
     EE_C_future = pi_E_C*CEE_C';
     EE_U_future = (EE_F_future*(pop_F/(pop_F+pop_C))+EE_C_future*(pop_C/(pop_F+pop_C)));
     
     pi_E_comp(1)=pi_E_R(1)*pop_R+pi_E_F(1)*pop_F+pi_E_C(1)*pop_C;
     pi_E_comp(2)=pi_E_R(2)*pop_R+pi_E_F(2)*pop_F+pi_E_C(2)*pop_C;
     pi_E_comp(3)=pi_E_R(3)*pop_R+pi_E_F(3)*pop_F+pi_E_C(3)*pop_C;
     pi_E_comp(4)=pi_E_R(4)*pop_R+pi_E_F(4)*pop_F+pi_E_C(4)*pop_C;
     pi_E_comp(5)=pi_E_R(5)*pop_R+pi_E_F(5)*pop_F+pi_E_C(5)*pop_C;
     
     pi_E_R_future(1) = pi_E_R*Pr_ee_R(:,1);
     pi_E_R_future(2) = pi_E_R*Pr_ee_R(:,2);
     pi_E_R_future(3) = pi_E_R*Pr_ee_R(:,3);
     pi_E_R_future(4) = pi_E_R*Pr_ee_R(:,4);
     pi_E_R_future(5) = pi_E_R*Pr_ee_R(:,5);
     
     pi_E_F_future(1) = pi_E_F*Pr_ee_F(:,1);
     pi_E_F_future(2) = pi_E_F*Pr_ee_F(:,2);
     pi_E_F_future(3) = pi_E_F*Pr_ee_F(:,3);
     pi_E_F_future(4) = pi_E_F*Pr_ee_F(:,4);
     pi_E_F_future(5) = pi_E_F*Pr_ee_F(:,5);    
     
     pi_E_C_future(1) = pi_E_C*Pr_ee_C(:,1);
     pi_E_C_future(2) = pi_E_C*Pr_ee_C(:,2);
     pi_E_C_future(3) = pi_E_C*Pr_ee_C(:,3);
     pi_E_C_future(4) = pi_E_C*Pr_ee_C(:,4);
     pi_E_C_future(5) = pi_E_C*Pr_ee_C(:,5); 
     
     pi_E_future(1) = pi_E_R_future(1)*pop_R+pi_E_F_future(1)*pop_F+pi_E_C_future(1)*pop_C;
     pi_E_future(2) = pi_E_R_future(2)*pop_R+pi_E_F_future(2)*pop_F+pi_E_C_future(2)*pop_C;
     pi_E_future(3) = pi_E_R_future(3)*pop_R+pi_E_F_future(3)*pop_F+pi_E_C_future(3)*pop_C;
     pi_E_future(4) = pi_E_R_future(4)*pop_R+pi_E_F_future(4)*pop_F+pi_E_C_future(4)*pop_C;
     pi_E_future(5) = pi_E_R_future(5)*pop_R+pi_E_F_future(5)*pop_F+pi_E_C_future(5)*pop_C;
     pi_E_future = pi_E_future./sum(pi_E_future);
     
% Consumption:
     C_M=(QM-xiH_F*pop_F-xiH_C*pop_C);
     C=QA+p_M*C_M+p_SR*QS_R+p_SU*QSU;
     C_A_share=QA/C;
     C_S_share=(p_SR*QS_R+p_SU*QSU)/C;

% Population distribution:
     sum_pop = pop_R+pop_F+pop_C;
     pop_R = pop_R/sum_pop;
     pop_F = pop_F/sum_pop;
     pop_C = pop_C/sum_pop;

     equilibrium = transpose([pop_R pop_F pop_C l_A l_M l_S...
                        pi_E_future(1) pi_E_future(2) pi_E_future(3) pi_E_future(4) pi_E_future(5)]);
     equilibrium_table = transpose([p_SR p_SU p_M wu wrR wrU wnU  0 xiH_F*p_M xiH_C*p_M ...
                QA_eq QS_R QSU QM GDP VA_A VA_SR VA_SU VA_SR+VA_SU VA_M GNI PYS_u PYS_rR PYS_rU PYS_nU...
                l_SR l_SU l_S l_M l_A pop_R pop_F pop_C...
                EE_R,EE_F,EE_C EE_U ...
                EE_R_future,EE_F_future,EE_C_future, EE_U_future ...
                pi_E_R(1) pi_E_R(2) pi_E_R(3) pi_E_R(4) pi_E_R(5) ...
                pi_E_F(1) pi_E_F(2) pi_E_F(3) pi_E_F(4) pi_E_F(5) ...
                pi_E_C(1) pi_E_C(2) pi_E_C(3) pi_E_C(4) pi_E_C(5) ...
                pi_E_U(1) pi_E_U(2) pi_E_U(3) pi_E_U(4) pi_E_U(5) ...
                pi_E_future(1) pi_E_future(2) pi_E_future(3) pi_E_future(4) pi_E_future(5)]);